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RELATION  OF  THE  ONE-PHASE  STEFAN  PROBLEM 
TO  THE  SEEPAGE  OF  LIQUIDS  AND 
ELECTROCHEMICAL  MACHINING 


Joel  C.  W.  Rogers 
1.  Formulation  of  the  Problems 

For  several  reasons,  I  have  recently  been  Interested  In  the 
problem  of  seepage  of  liquids  through  a  dam.  First,  my  experi¬ 
ence  with  algorithms  for  computing  generalized  solutions  of  other 
hydrodynamic  free  boundary  problems  (water  waves.  Ref.  9)  has  nat¬ 
urally  led  me  to  the  question  as  to  the  applicability  and  useful¬ 
ness  of  such  methods  for  the  dam  problem.  Second,  the  "Balocchl 
transformation"  which  has  been  used  effectively  to  solve  some  free 
boundary  problems  of  seeping  liquids  is  reminiscent  of  another 
transformation,  relating  the  steady  state  of  a  one-phase  Stefan 
problem  to  a  time-dependent  free  boundary  problem  arising  In  the 
theory  of  anodic  smoothing  (Ref.  8),  and  the  possibility  that  the 
two  transformations  might  have  a  common  origin  has  Intrigued  me. 

A  third  reason  Is  that.  In  the  second  announcement  of  this  Inten¬ 
sive  seminar,  I  was  billed  as  discussing  "numerical  methods"  and 
the  dam  problem,  which  Is  related  through  the  Balocchl  transforma¬ 
tion  to  other  free  boundary  problems  In  whose  numerical  solution  I 
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have  participated  (Refs.  2-5),  seemed  like  an  apt  problem  to 
discuss.  Finally,  there  Is  no  better  forum  In  which  1  can  discuss 
this  subject  than  this  one  In  Pavla,  where  so  much  work  on  the  dam 
problem  has  been  done. 

The  format  of  this  paper  Is  as  follows:  In  the  remainder  of 
this  Introductory  section  we  will  formulate.  In  a  cursory  fash¬ 
ion,  the  free  boundary  problem  associated  with  the  seepage  of  a 
liquid,  as  well  as  variants  of  the  one-phase  Stefan  problem;  for 
these  we  will  also  provide  algorithms  which  have  been  used  to 
effect  their  numerical  solution.  The  next  section  will  contain 
a  discussion  of  the  time-dependent  porous  flow  problem.  Its  rela¬ 
tion  to  the  time-dependent  anodic  smoothing  problem,  and  some 
comments  about  the  numerical  solution  of  the  seepage  problem.  In 
a  third  section  we  examine  the  solution  of  the  steady  state 
porous  flow  problem  for  rather  general  dam  shapes  In  greater 
detail;  we  propose  an  algorithm  for  the  solution  of  the  problem 
and  make  some  tentative  remarks  pertinent  to  the  question  of  error 
estimates  for  the  approximate  solution  generated  by  the  algorithm. 
A  final  section  Indicates  briefly  a  possible  generalization  of 
this  work.  Owing  to  time  limitations,  the  results  I  report  In 
this  paper  are  not  as  complete  as  I  would  like  them  to  be,  and  I 
apologize  In  advance  for  this  Incompleteness.  Nevertheless,  If 
the  perspective  offered  In  this  paper  usefully  complements  other 
work  on  the  dam  problem  and  thereby  contributes  to  Its  solution, 

I  will  consider  the  paper  to  have  served  Its  purpose. 
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By  the  one-phase  Stefan  problem  In  a  bounded  region  (Refs. 
5,  8),  we  Mill  mean  the  problem  of  finding  the  solution  of  the 
equation 

ut  «  Af(u)  ,  x  e  0  ,  0  <  t  <  »  ,  (1.1a) 

subject  to  boundary  conditions 

Lf(u)  “0  ,  x  e  30  ,  0  <  t  <  »  ,  0.1b) 

and  Initial  conditions 

u(x,0)  »  u0  ,  x  £  D  (1  .lc) 

where  L  Is  an  appropriate  linear  operator  and  the  function  f  Is 
given  by 

f(u)  =  max(u  -1,0)  .  (1 .2) 

Under  most  conditions  of  interest,  the  following  algorithm 
suffices  for  the  calculation  of  u(x,t)  (Refs.  4,  5,  8): 

u(n-r)  ~  un  ,  (1.3a) 

u°  -  u0  ,  (1.3b) 

un+1  -  un  +  (Si(x)  -  l)f(un)  ,  (1.3c) 

where  the  operator  Si(x)  Is  defined  by 

C(t)  s  Si(t)?o  (1.4a) 

when  C  satisfies 

,  x  c  0  ,  0  <  t  <  t  (1 .4b) 

U  »  0  ,  x  £  30  ,  0  <  t  <  t  ,  (1 ,4c) 

5(0)  =  50  .  x  £  0  .  (1 ,4d) 

We  shall  refer  to  the  algorithm  (1.3)  as  algorithm  I,  A  steady 
state  u  exists,  and  algorithm  I  may  be  used  to  calculate  IT. 
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A  variant  on  the  one-phase  Stefan  problem  which  has  a  more 
direct  connection  to  the  dam  problem  Is  the  following  set  of  con¬ 
ditions  satisfied  by  the  function 


v(x,t)  ■  |  f(u(x,t')dt' 

(1.5) 

when 

Uo (x )  >_  1  for  x  e  fi(0)  c  d  , 

(1.6a) 

uq (x )  =  0  for  x  £  D  -  n(o)  : 

(1.6b) 

v^  «  Av  +  Uo  -  1  ,  x  e  fl(t)  ,  0  <  t  <  « 

,  (1.7a) 

,|3D(t)-30  '  0  '  0  4  1  ‘  *  • 

(1.7b) 

vw<3n(t)-ao  '  0  •  0  <  1  ‘  "  • 

(1.7c) 

Lv  a  0  ,  X  e  3  D  ,  0<t<»  , 

(1.7d) 

v(x,0)  ■  0  ,  X  e  D  , 

(1.7e) 

which 

has  the  property  that 

supp  v(x,t)  *  n(t)  c  D  , 

(1.8a) 

n(ti )  a  n(t2)  for  ti  >_  t2 

(1.8b) 

v  and  Vv  are  continuous  In  D 

(1.8c) 

Obviously,  we  may  use  (1.5)  and  algorithm  I  to  calculate 

v .  How- 

ever. 

problem  (1.7)  Is  of  a  type  which  occurs  In  the  theory  of 

oxygen  transport  In  tissue,  and  the  following  algorithm  may  be 

used  to  calculate  v(x,t)  (Ref.  3): 

v(nr)  ~  vn  , 

(1.9a) 

v°  •  0  , 

(1.9b) 

vn+1  ■  max(S2(T)vn  -  t,0)  , 

(1.9c) 

where 

S2(t)  Is  defined  by 

C(t)  s  S2(t)Cq 

(1.10a) 

4 


when  t  satisfies 


*  A£  +  U0  ,  X  £  D  ,  0  <  t  <  T  , 

*  (1.10b) 

U  *  0  ,  x  «  3 D  ,  0  <  t  <  t  ,  (1 ,10c) 

5(0)  «  50  .  x  e  D  .  (l.lOd) 

We  refer  to  the  algorithm  (1.9)  as  algorithm  II.  For  this  algo¬ 
rithm  an  explicit  error  estimate  Is  available  (Ref.  3).  It  Is 

XS“P0  | v(x,nt)  -  vn(x)|  m  .  (1.11) 


In  a  system  of  coordinates  (x,t)  In  which  the  gravitational 
acceleration  Is  unity,  an  Irrotatlonal  and  Incompressible  flow  In 
a  porous  medium,  assumed  to  be  an  open  set  D%  Is  given  by  a  veloc¬ 


ity  potential  ?  satisfying  the  equations 

?ya-Z-p-a?-  4(7?) 2  ,  x  e  n(t)  c  D  ,  0  <  t  <  ®  , 

1  (1.12a) 

A?  =  0  ,  X  e  U(t)  ,  0  <  t  <  “  ,  (1.12b) 

p  =  0  In  D  -  n(?)  ,  (1 .12c) 

where  p  Is  the  pressure  and  a?  Is  the  potential  for  a  frictional 

drag  force  on  the  flow.  In  the  dam  problem  In  N  dimensions, 

D  c  RN"1x(0,»)  ,  (1.13a) 

and  with  3 D  decomposed  as 

3 D  =  3D0  u  3Dp  u  3 Dj  ,  (1.13b) 

the  following  boundary  conditions  are  Imposed: 

p  ■  0  on  3Z>o  »  (1 .13c) 

p  ■  zR  -  z  »  hydrostatic  pressure, 

or  more  generally  p  ■  pR  >  0  on  3 z?R,  (1.13d) 

n  •  7?  Is  prescribed  on  3Z? j  .  (1.13e) 


Making  the  further  decomposition  3£>j  ■  30^  u  3£)j2  where 

S0I!'  su«>  ,  we  require  that 

3  D  « 

1  3fl(t)  oZDr  u  3Dj2  .  (1.13f) 

The  only  case  of  physical  Importance  Is  the  case 

p  >_  0  throughout  D  .  (1.1 3g ) 

This  constraint  should  be  viewed  as  a  restriction  on  the  types  of 
Inflows  and  outflows  which  may  be  prescribed  on  3 Dj  In  (1.1 3e ). 
A  case  frequently  considered  Is  that  for  which  3 Is  Imperme¬ 
able:  3 z?l2  »  0.  In  the  literature,  many  authors  have  made  the 
further  assumption  that  3 a  B0,  where 

B  -  ( (x ,z )  €  Ci(D)|z  =  z0)  .  (1.14) 

3 Z?R  Is  that  portion  of  3 D  In  contact  with  a  reservoir  of  pre¬ 
scribed  surface  height  zR,  and  z  <  zR  on  3DR.  We  assume  that  p, 
given  on  3 D0  u  3DR  by  (1.13c)  and  (1.13d),  has  bounded  derivatives 
on  3D0  u  3Dr.  The  boundary  of  n(t)  Is  decomposed  as 

3£5  =  3fl  n  W  u  3flf  .  (1  .15) 

Boundary  conditions  on  30^  are 

p  =  0  ,  x  e  3n^(t)  ,  0  <  t  <  «  ,  (1.16) 

and  points  x(t)  on  3n^(t)  move  with  the  velocity 


An  Initial  value  problem  for  the  flow  in  D  would  have  as  Its 
object  the  determination  of  a  solution  ($(t),n(t))  of  (1.12)- 
(1.17)  from  the  Initial  data  (l>(0)  ,n(0)) .  However,  this  problem 
will  not  generally  possess  a  unique  solution.  The  reason  Is  that. 
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for  an  Irrotatlonal  Incompressible  flow  In  a  region  0,  one  Is 

free  to  specify  not  only  the  velocity  potential  where  the  density 

Is  nonzero,  but  also  the  density  where  the  flow  enters  0.  In  our 

problem  the  velocity  field  Is  given  by  $  where  It  Is  specified, 

and  the  density  (which  essentially  takes  on  only  the  values  0  and 

1  )  is  given  by  specifying  the  region  n(t)  occupied  by  fluid. 

Accordingly,  given  a  flow  region  n(t),  we  are  free  to  allow  fluid 

to  enter  anywhere  on  3£0-  3n(t),  thereby  augmenting  the  flow 

region.  The  restriction  (1.12c)  only  limits  the  types  of  such 

Inflows  to  masses  of  fluid  In  free  fall  In  a  frictional  medium. 

Such  a  flow  can  be  found  by  solving  (1.12a)  with  p  =  0,  yielding 

^(x.t)  *  e-«(t-t°)Y  .  x  _  i/a(i  .  e“a(t-t0)jz  +  constant 
~  ~  (1.18a) 

In  a  region  n0(t)  moving  with  the  fluid  with  the  spatially 

constant  velocity 

•f 

v?0  *  Ye-<x(*"*o)  -  ^(1  -  e"a(*'*o))  ,  (1.18b) 

a 

and  making  its  appearance  at  time  t0  In  d  at  points  of 

3£>o  -  3n(t0)  for  which  n  •  y  <  0.  Denote  the  augmented  region  by 

n+(t)  =  !!(!)  u  n0(!)  .  (1.19) 

It  follows  from  our  discussion  that  the  velocity  potential  $+Ct) 

~+/~,  **(f0) 

defined  In  n  (t)  will  have  the  property  that  — — - <  0  for  some 

point  on  3n+(?0)  n  3 0O.  We  can  make  the  solution  of  (1.12)- 
(1.17),  subject  to  given  Initial  conditions,  unique,  by  requiring, 
for  a  given  decomposition 

30O  *  301  u  302  (1.20) 


of  30o «  that  for  X  >  0 
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(1.21a) 


r 


1  0  ,  X  e  3ft(t)  n  3D*  0  .21a) 

and 

3fi(t)  3  3 D2  .  (1 .21b) 

In  flow  through  a  porous  medium,  the  case  of  Interest  Is 
that  In  which  the  frictional  drag  coefficient  a  Is  so  large  that 
the  dependent  variables  describing  the  flow  change  Insignificantly 
over  the  time  1/a.  We  may  think  of  1/a  as  a  "relaxation  time" 
during  which  the  flow  assumes  the  asymptotic  form  associated  with 
the  limit  a  -*•  «.  The  actual  asymptotic  parameter  Is  a2a,  where 
a  Is  a  characteristic  dimension  of  the  dam  D. 

Because  of  the  large  drag  force,  It  Is  clear  that  we  will  be 
dealing  with  a  very  slow  flow,  and  that  In  the  limit  as  a  » 
there  will  have  to  be  some  balance  In  (1.12a)  between  the  drag 
potential  a$  and  the  potential  of  the  gravitational  and  pressure 
forces,  z  +  p.  Thus,  we  shall  work  with  the  dependent  variable 

♦  -  a?  .  (1 .22) 

Since  the  relaxation  time  Is  0(l/a),  we  may  Investigate  the 

phenomenon  of  relaxation  by  working  with  the  Independent  variable 

t  *  at  (1 .23) 

Then  at  a  point  on  3^,  using  (1.16),  (1.17),  and  (1.12a),  we  get 
for  the  rate  of  change  with  time  of  $  for  a  point  moving  with  the 
fluid, 

5S-1-*  0-24) 

at 

while  from  (1 .17) 
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0.25) 


$  a  -  z  ,  x  e  anf(t)  ,  0  <  t  <  ».  (1.30) 

As  a  +  «,  we  can  use  (1 .12a),  (1 .12b)  and  (1 .22)  to  get 


P  3  -  4>  -  z  ,  x  e  n(t) 

0  <  t  <  “  , 

(1.31a) 

AP  =  0  ,  x  £  n(t)  , 

0  <  t  <  »>  , 

(1.31b) 

p  »  0  on  3D0  , 

(1.31c) 

p  =  pR  >  0  on  3Cr 

» 

(1 . 31  d ) 

n  •  vp  prescribed  on  3Cj  n 

30(t)  , 

(1.31e) 

p  =  0  on  3flf(t) 

» 

(1 .31  f) 
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(1.31g) 


f 

and  3 (if  moves  with  the  velocity 

-  Vp  -  k 

As  before,  we  restrict  our  attention  to  problems  for  which 

p  >_  0  ,  x  c  D  ,  0  <  t  <  00  ,  ( 1  . 31  h ) 

and  we  require 

3fi(t)  =  3/>r  u  9D  .  (1  .311) 

Finally,  to  Insure  uniqueness,  we  use  the  conditions  (1.21)  in 
the  form 

-  >  k  •  n  ,  x  c  an(t)  n  jd,  ,  0  <  t  <  «, 

9n"  (1 .31 j) 

3n(t)  3  3 d2  »  0  <  t  <  oo  .  (1.31k) 

For  the  rest  of  this  paper,  we  will  consider  (1.31)  to  be  the 

time-dependent  version  of  the  dam  problem.  This  free  boundary 

problem  Is  very  similar  to  the  problem  for  the  electrostatic 

potential  In  the  theory  of  anodic  smoothing  (Ref.  8),  and  In  fact 

the  two  problems  would  be  Identical,  were  It  not  for  the  term  -k 

In  ( 1 . 31 g ) ,  and  the  different  boundary  conditions  Imposed  for  p 

on  the  "cathode"  surface  3£>R  u  3£?j. 

Remark  1.1;  If  we  consider  the  evolution  of  a  front  In  an  Ini¬ 
tially  dry  porous  medium  abutting  reservoirs.  It  will  follow  from 
(1.31g)  that  the  Initial  velocity  of  the  front  is  unbounded.  This 
reflects  an  error  in  the  approximation  (1.30),  which  only  holds 
after  a  time  t  *  0(1 /a),  and  also  an  error  In  using  the  hydro¬ 
static  assumption  (1.31d)  for  the  reservoirs,  as  this  breaks  down 
during  the  Initial  phase  of  the  flow. 


In  the  beginning  of  the  Introduction  to  this  paper,  we 
mentioned  the  possible  applicability  of  algorithms  for  other 
hydrodynamic  free  boundary  problems  to  the  solution  of  the  time- 
dependent  dam  problem.  Essentially  In  those  more  general  algo¬ 
rithms  the  free  boundary,  the  motion  of  which  Is  given  by  (1.29), 
Is  determined  by  looking  for  contours  p  =  constant  where  p 
satisfies  an  equation  like 

Pj  +  5^  1  Vp  a  0 

and  p  has  a  jump  at  the  free  boundary.  The  more  general  hydro- 
dynamic  algorithms  do  not  appear  to  have  any  utility  for  the 
solution  of  the  dam  problem. 

2.  Comments  on  the  Time-Dependent  Dam  Problem 

A  monotonicity  result  can  be  deduced  Immediately  for 
solutions  of  the  dam  problem  (1.31).  We  formalize  It  as  a  lemma. 

Lemma  2.1:  Suppose  we  have  two  solutions  pi  and  p2  of 
(1.31)  for  which  PiI3DrU3D()  1  P2 1 3£?Ru3Z?0  ’  Sn131  UzJj  -an^UzJj’ 

and  «x(0)  o  n2(0).  Then  fii(t)  =>  n2(t). 

Proof:  If  for  some  t,  we  have  nj(t)  n2(t),  then  at  some 
ti  <  t  a  point  of  an^(ti)  first  passed  through  3ft^(ti).  Then 
ni(ti)  a  02(ti)  and  3  a  point  P  e  an^tj)  n  3n^(t1)  for  which 


11 


‘  a7T(p»tl)  <  *  37T(p*t»)  • 

But  a(pi  -  p2)  ■  0  In  n2(ti),  Pi  -  P2  1  0  on 

an2(ti)  n  (3  D0  u  3£>R),  ~<Pi  -  P2 )  >  on  3ft2(tx)  n  3Z)j,  pj  -  p2  >  0 

on  3ft.  (tj),  and  (pj  -  p2)(P,tj)  ■  0.  Hence  we  get 
“2 


3P1  3  P2 

“ar{p,tl)-"  3H_{p'tl)  ’ 

giving  a  contradiction  and  proving  the  lemma. 

It  follows  right  away  that  If  the  prescribed  boundary  data 
are  time-dependent  and  If  ft(t)  Is  Increasing  (decreasing)  at  any 
time  t0,  that  Is,  n(t0  +  0  3  ft(t0)(ft(t0  +  0  c  ft(t0))  for  some 
6  >  0,  then  ft(t  +  6)  =  ft(t)(ft(t  +  <$)<=  ft(t))  for  all  t  >_  t0 .  By 
picking  regions  ft(0)  for  which  one  can  determine  a  priori  that 
-  |£  -  £  «  n^0(^0)  everywhere  on  3ft^(0),  one  can  then  estab¬ 
lish  the  existence  of  a  steady  state  solution,  assuming  that  a 
solution  to  the  time-dependent  problem  exists. 

The  following  stability  result  Is  also  easily  deduced. 

Lemma  2.2:  Suppose  we  have  two  solutions  pi  and  p2  of  (1.31) 
with  the  same  boundary  data  prescribed  on  3 D.  Then 

gtl(Mt)  -  fl2(t))  u  (n2(t)  -  ft!(t))|  <  0  .  (2.1) 


Proof:  Without  loss  of  generality,  we  can  have 
fti(0)  3  ft2(0),  and  all  we  need  to  verify  Is  that 


atlMt)  -  ft2(t)i 


t-o  -  From  ^ 


12 


jtlMO  -  Mt)|  t-0 


f  -  (3^-  -  k  .  n)dS 

J»Qf  (0)  V3n  / 

f  /3  P2  ♦  \ 

+  Li^for3"”  *  k  *  n'ds 

f  (stt  +  *  *  n)ds 

J  30i  (0  )n3£)' 3n 

[  /  3  P2  -*-  \ 

(rr-  +  k  •  n  dS 

■’3fl2(0)n3Zr 90  ' 


(2.2) 


Since  3£^1  (0)  n  3  Do  an2(0)  n  3D  3  supp  Pflgp.  Pi  t.  P2  L  0 
Vx  e  D,  and  g-^Pi  -  p2)  *  0  or  pi  -  p2  *  0  on  3R2(0)  n  3 z>,  it 
follows  from  (1.311)  and  (1 .31 j)  that 


3tlni(t)  -  n2(t)| 

♦J 


t-o 


f 


dSl2(o)ndD 
(*  Pi 


3  ( Pi  "  P2) 
3rT 


-dS 


(3^1 (o)-3fi2(o))n3Dj 


/  »Pl  -  \ 

\3n“  +  k  *  n  ds/l  0  •  (2. 


3) 


(2.3)  may  be  combined  with  lemma  2.1  to  prove  (2.1). 

Remark  2.1:  It  seems  likely  that.  If  we  restrict  ourselves 
to  components  of  fii  and  n2  which  are  connected  to  3 o0  u  3 z?R,  we 
can  get  a  stronger  relation  than  (2.3),  and  prove  decay  of 
|Ri (t)  -  n2(t)|  to  0  as  t  +  «.  Indeed,  In  these  circumstances 
It  appears  that  for  some  constant  b  >  0  we  can  show  that 
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3tlni(t)  -  n2(t)|  <  -  b|nt(t)  -  n2(t)|  , 

In  which  case 

IMt)  -  n2(t)(  <_  e“bt|ai(0)  -  n2(0)|  .  (2.4) 

I  have  been  Informed,  In  the  course  of  writing  this  up,  that  a 
counterexample  has  been  given  by  Alt  for  components  of  Jli  and  n2 
which  are  not  connected  to  3 D0  u  a (Ref.  1). 

An  Immediate  consequence  of  lemma  2.2  Is  the  following:  If 
the  prescribed  boundary  data  are  time-independent,  the  Initial 
rate  of  change  of  n(t)  Is  the  largest,  that  Is, 

|(n(t  +  5)  -  fi(t))  u  (n(t)  -  n(t  +  6 ) ) |  i  |(n(s)  -  n(o)) 

u  (0(0)  -  0(6))l  . 

(2.5) 

Returning  to  direct  consideration  of  the  problem  (1.31),  we 
recall  our  remark  regarding  Its  similarity  to  the  anodic  smooth¬ 
ing  problem.  We  found  (Ref.  8)  that  the  solution  of  the  latter 
problem,  at  a  time  t,  could  be  found  directly  by  finding  the 
steady  state  of  problem  (1.7)  with  u0  now  dependent  on  the  param¬ 
eter  t.  Making  some  alterations  and  specifications  In  the  steady 
state  version  of  (1.7)  In  order  that  the  result  will  conform  to 
the  sort  of  boundary  value  problem  we  have  been  considering  here, 
we  note  that  the  solution  to  the  "anodic  smoothing"  problem  at 
time  t  can  be  found  from  the  solution  of  the  following  elliptic 
free  boundary  problem: 
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AC  - 

'  0  , 

x  e  n(0)  c  d  ,  0  <  t  <  «*> 

,  (2.6a) 

AC  -  1 

,  X  6 

n(t)  -  n(0)  c  d  ,  0  <  t  < 

*  » 

(2.6b) 

C  >  0 

,  X  £  D  ,  0  <  t  <  «  , 

(2.6c) 

C 

-  0  , 

X  e  3£>0  ,  0  <  t  <  « 

( 2 . 6d ) 

C  - 

tpR  >  0 

,  x  e  3  Z?R  ,  0  <  t  <  * 

,  (2.6e) 

3C 

3n 

.  t^ 

3n 

,  xeSDj,  0<t<“  , 

(2.6f) 

C  5 

0  ,  X  e-  D  -  n(t)  , 

(2.6g) 

C  and 

vc  are  continuous  In  D 

(2.6h) 

Differentiation  of  each  of  the  conditions  (2.6)  with  respect 
to  t  shows  that  ct  satisfies  the  problem  (1.31),  except  for  the 
gravitational  term  -k  In  (1 .31 g ) .  If  we  transform  coordinates 
and  view  the  problem  from  a  system  which  Is  moving  downward  with 
unit  velocity,  the  term  -k  will  disappear,  but  then  D  will  appear 
to  be  moving  upward  with  unit  velocity.  Thus,  the  problem  is 
that  of  anodic  smoothing  when  the  "cathode"  Is  In  uniform  motion. 

■f 

It  should  be  noted  that  the  term  -k  In  (1 .31 g )  effects  a  uni¬ 
form  translation  of  the  free  surface  and  has  no  effect  on  the 
character  of  Int(3fl^,)  with  regard  to  Its  differential  geometric 
properties.  Thus,  the  observed  relation  of  the  solution  of  the 
time-dependent  dam  problem  to  the  solution  of  a  sequence  of  steady 
state  Stefan  problems  In  a  translating  region  may  be  useful  pri¬ 
marily  for  what  can  then  be  deduced  about  the  regularity  of  the 
free  surface  from  known  results  for  the  steady  state  Stefan 
problem,  and  not  for  practical  calculation. 
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Nevertheless,  we  can  exploit  this  relation  to  suggest  a 
numerical  approach  to  the  time-dependent  dam  problem:  Viewed 
from  a  system  translating  uniformly  downward  with  unit  speed, 
each  point  P  e  R  fixed  In  space  Is  moving  upward  uniformly  with 
unit  speed,  and  moves  through  the  distance  t  In  the  time  interval 
[nt,(n  +  1 ) t] .  Let  us  denote  by  cn(r)  the  points  thus  swept  out 
by  a  set  r  c  RN.  We  define 

Dn  i  Sn(D)  ,  (2.7) 


a  region  which,  for  each  n.  Is  fixed  with  respect  to  the  moving 
system.  Given  n(0)  c  d ,  we  construct  fin  according  to  the  follow¬ 
ing  prescription.  First  we  define 


h°  =  o°(n(0))  c  5 0  . 

(2.8) 

Then,  given  n",  n  0,  we  find  nn+1  by  solving 

AC  =  0  ,  x  e  fin  c  Dn 

(2.9a) 

ac  *  l  ,  x  e  n"+1  -  nn  c  bn  , 

(2.9b) 

C  _>  0  ,  X  c  Dn  , 

(2.9c) 

C  a  0  ,  X  e  on(&00)  n  dDn  , 

(2.9d) 

C  *  tpR  >  0  ,  X  e  on(3Z?R)  n  3 5°  , 

(2.9e) 

.  X  «;"(»£>,)  n  30"  . 

(2.9f ) 

e  i  o  ,  X  «  tr  -  s?+1  . 

(2.9g) 

C  and  ?c  are  continuous  In  Dn 

(2.9h) 

nn+l  Is  determined  from 

6ntl  3  n"*1  n  d"+1  . 

(2.10) 

Finally,  the  set  of  points  fin  c  5n,  when  viewed  from  the  fixed 

*n 

system  at  time  nr.  Is  a  set  n  ,  and  we  define 
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nn  H  n*n  n  D  .  (2.11) 

Viewed  from  the  system  at  rest,  the  region  Dn  Is  seen.  In 
the  time  Interval  [nx,(n  +  1)t),  as  a  time-dependent  region  D*{ t) 
Independent  of  n  and  given  as  follows: 

D*{ t)  “  a*(t ){D)  (2.12a) 

where  for  any  set  r  c  R  ,  o  (t)(r)  Is  defined,  for 
t  c  [nx,(n  +  1 )t),  by 

a  (t)(r)  =  | (x,z)| (x,z' )  c  r  for  some 

z'  £  [z  +  t  -  (n  +  1)t,z  +  t  -  nt]| 

*n  ,  (2.12b) 

n  will  then  be  the  region  fl(nx)  solving  the  problem  (1.31), 


with  the  following  modifications: 

n(0)  Is  replaced  by  o*(0)n(0)  ,  (2.13a) 

d  Is  replaced  by  D*[t)  ,  (2.13b) 

3D0  1s  replaced  by  o*(t)(az?0)  n  3D*(t)  ,  (2.13c) 

3Dr  Is  replaced  by  o*(t)(3£>R)  n  az?*(t )  ,  (2.13d) 

3Dj  Is  replaced  by  a*(t)(3Z?I)  n  3D*(t)  .  (2.13e) 


If  3 d  Is  at  all  smooth  (having,  say,  piecewise  bounded 
curvature),  the  same  sort  of  smoothness  (up  to  boundedness  of 
curvature)  will  characterize  3o*(t).  Accordingly,  a  priori  regu¬ 
larity  results  will  be  derivable  for  p  satisfying  the  modified 
version  of  (1.31),  and  since  d  (t)  :>  ov t,  one  may  deduce  the 
deviation  between  the  actual  values  attained  by  p  and  on  3P 
for  the  solution  of  the  modified  problem,  and  the  values  required 
In  the  original  problem  (1.31).  To  obtain  error  estimates  for 
the  solution  to  the  modified  problem,  as  opposed  to  the  solution 
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of  (1.31),  one  need  only  determine  a  domain  n"  known  a  priori 
to  be  contained  for  all  time  In  the  union  of  the  supports  of  p 
for  the  modified  and  original  problems,  a  determination  easily 
made  with  the  help  of  lemma  2.1,  then  find  the  capacities  rela¬ 
tive  to  n”  of  various  subsets  of  a D0  u  and  multiply  them  by 

the  corresponding  errors  In  the  values  of  p  found  for  those  sub¬ 
sets,  Integrate  the  error  In  over  a Dj,  and  finally  make  use  of 
lemma  2.2.  For  cases  of  Interest,  we  anticipate  an  error,  as 
measured  by  the  volume  contained  between  the  two  determinations 
of  the  free  boundary,  proportional  to  t  for  finite  t.  If  the 
expectation  of  remark  2.1  Is  borne  out,  such  an  error  estimate 
will.  In  fact,  hold  uniformly  for  all  time. 

We  recall  that.  If  p  satisfies  a  problem  like  (1.31)  for  p, 
except  for  the  absence  of  the  term  -k  In  (1.31g),  then 

e(x,z,t)  -  |  p(x,z,t')dt'  (2.14) 

satisfies  the  elliptic  free  boundary  problem  (2.6).  This  Is  the 
transformation  used  to  solve  the  anodic  smoothing  problem 
directly.  (It  Is  somewhat  reminiscent  of  the  transformation  (1.5) 
relating  variants  of  the  Stefan  problem.)  Since  the  term  -k  in 
(1.31g)  has  been  seen  to  vanish  In  a  coordinate  system  moving 
downward  with  unit  velocity,  the  natural  adaptation  of  (2.14)  to 
the  dam  problem  Is  the  transformation 

u(x,z,t)  ■  p(x,z  +  t  -  t',t')dt'  .  (2.15) 

When  a  steady  state  p(x,z)  has  been  achieved  as  t  +  ®,  (2.15) 
becomes 
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(2.16) 


u(x,z)  ■  j  p(x,z  +  t')dt* 

Jo 

which  Is  precisely  the  Balocchl  transformation  (Ref.  2). 

Note  added  In  proof:  I  have  learned  since  doing  this  work 
that  (2.15)  Is  exactly  the  same  as  the  transformation  Introduced 
by  Torelll  In  a  paper  by  Friedman  and  Torelll  (Ref.  7). 

The  formal  transformation  (2.15)  Is  not  defined  when 
(x,z+t-t')^0  for  some  t'  e  [0,t].  (The  extension  of  p  to 
D  -  n(t)  through  the  convention  that  p  =  0  there  Is  quite  natu¬ 
ral.)  Although  a  version  of  (2.15)  can  be  given  which  Is  satis¬ 
factory  for  the  time-dependent  case  In  a  fairly  general  region  D , 
we  shall  not  consider  It  here.  In  the  next  section,  for  the 
steady  state  case,  a  version  of  (2.16)  appropriate  for  more  gen¬ 
eral  regions  D  will  be  considered.  For  the  present,  let  us 
restrict  ourselves  to  the  case  In  which  D  Is  the  cylinder 
(cf.  (1.14)) 

D  »  B0x(0,»)  ,  (2.17) 

and  also 


3Z?r  -  B„  ,  (2.18a) 

3Z>0  u  3 Dr  »  3B0x(0,»)  .  (2.18b) 

In  this  case  (x,z  +  t  -  t')  c  Z?vt'  e  [0,t].  The  time- 
dependent  dam  problem  (1.31)  Is  then  equivalent  to  the  following 
problem  for  u: 


u  and  ?u  are  continuous  In  v  , 
u  s  0  ,  x  €  D  -  n(t)  ,  0  <  t  <  • 


(2.19a) 

(2.19b) 
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u  >_  0  ,  x  e  E>  ,  0  <  t  <  «  ,  (2.19c) 

u  given  ,  x  e  3B0x(0,«>)  ,  0  <  t  <  «  ,  (2.19d) 

ut  -  u  given  ,  x  £  B0  ,  z  ■  0  ,  0  <  t  <  <*>  , 

tZ  ZZ  (2.1 9e) 

(1  If  (x,z)  e  fl(t)  and  (x,z  +  t)  J  n(O) 

Au(x,z,t)  *  |  (2.1 9f ) 

(o  otherwise 

When  a  steady  state  u(x,z)  Is  reached,  the  problem  (2.19)  reverts 
to  a  form  which  can  be  expressed  as  a  variational  Inequality  and 
which  has  been  considered  already  In  some  detail  (Ref.  2).  This 
problem  Is  very  close  to  the  steady  state  version  of  the  free 
boundary  problem  (1.7)-(1.8),  which  can  be  solved  by  either  algo¬ 
rithm  I  or  algorithm  II. 


3.  Steady  State 


Our  first  order  of  business  Is  to  give  a  meaningful  version 
of  the  transformation  (2.16)  for  a  fairly  general  region  D.  To 
this  end,  we  will  find  It  convenient  to  work  with  functions 
satisfying  homogeneous  conditions  on  3 D,  and  we  suppose  that  we 
can  write 


p  =  _  +  ii»  (3.1 ) 

where  ez  has  the  following  properties: 

supp(ej  -  C fc(  3Z?d)  -  C *.(  3 D.)  c  sf  ,  where  we  know  a  priori 

Z  K  12 

that  c  Int(n)  , 
(3.2a) 


A0Z  <_  0  V  x  c  D  , 

0Z  ■  0  ,  X  e  3 D0  , 


(3.2b) 

(3.2c) 


20 


.z  .  -  pR  <  0 

,  x  C  3£>R 

(3. 2d) 

rn*z  ■  •  In  '  • 

X  e  *Dl2 

(3.2e) 

T<?z  *  0  • 

X  e  30^ 

(3.2f) 

Obviously, 

,  from  (1 .31 )  and  (3.2) 

if*  ■  0  ,  x 

e  3O0  u  3 £>R  , 

(3.3a) 

• 

X  £  3Z?I2 

(3.3b) 

=  A0z  <_  0 

,  X  £  0  , 

(3.3c) 

♦  *  o  , 

X  £  30^.  , 

(3.3d) 

U*  -  *  • n  ■ 

X  £  30^  n  3ft 

(3.3e) 

As  noted  after  equation  (1.13g), 

we  have  i|>  >  0  on  30^  n 

30.  In 

addition. 

since  p  >_  0  In  D  and  p 

=  0  on  3R^,  we  get  from  (3.2) 

and  (1 .31 g )  that 

,  X  £  3ft^ 

(3.3f ) 

All  these 

conditions  Imply  that 

*  >  0  , 

x  £  ft 

(3.3g) 

It  Is 

now  natural  for  us  to 

make  the  extension 

♦  =  0  , 

x  £  d  -  n 

(3.4) 

Defining, 

for  all  x  e  RN"\ 

Z(x)  =  sup{  z  | 

(x,z)  £  D} 

(3.5) 

we  are  led,  tn  our  version  of  the  Balocchl  transformation,  to  con¬ 
sider  the  function 
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rZ(x) 

c(x,z)  -  ^(x,z*)dz'  .  (3.6) 

z 

We  will  rewrite  the  steady  state  problem  as  a  problem  for  c, 
In  a  form  for  which  a  solution  algorithm  readily  suggests  Itself. 
Before  we  proceed,  however.  It  will  be  useful  to  distinguish  sev¬ 
eral  cases  which  may  arise.  It  will  prove  helpful  for  us  to  make 
the  definitions,  for  x  £  RN_1 , 


(zft(x)>  ■  3{z|(x,z)  £  n)  , 

(3.7a) 

z  (x )  ■  sup(z  £  {zn(x)D 

(3.7b) 

(zf(x)l  «  (z| (x,z)  £ 

(3.8) 

Also  we  define  the  set 

{ZD(x)l  »  3{z[(x,z)  £  D } 

(3.9a) 

whose  members  are  written  as 

zi(x)  >  z2(x)  >  ...  >  zn(x) (x) 

(3.9b) 

If  Z(x)  =  «,  n(x)  Is  odd,  and  on  account  of  (3.3e)  elements 
of  {Zj,(x)}  can  only  lie  between  <*>  and  z^x),  or  between  z2l(x) 
and  ^2 1+1 (x )  for 

1  c  SD(x)  =  jl|1  <.1  <  ~  1  ,  (x,Z2l(x))  e  9DX  u  BCjjj 

(3.10a) 

If  Z(x)  <  «,  n(x)  Is  even,  and  elements  of  {z^(x)>  can  only  lie 
between  z  (x)  and  z2l{x)  for 

1  c  S^(x)  =  jl|l  <  1  <  ,  (x,z2l_1(x))  £  %D\  u  3 Z?  1 1  i  • 

(3.10b) 

If  n(x)  Is  0,  1,  or  2,  we  see  that  the  number  of  elements  of 
{zn(x)>  Is  0  or  2.  In  the  former  case  (z^(x)}  Is  empty;  In  the 
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latter  It  Is  empty  or  contains  at  most  one  element,  z*(x).  Note 
that  If  D  Is  convex,  we  have  that  n(x)  Is  0,  1,  or  2  for  all  x. 

We  shall  now  restrict  ourselves  to  the  case  where  for  each  x 
the  set  (Zj(x)}  contains  at  most  one  element,  z  (x),  and  also 
where  v x 

{(x,z)  e  Cl{D)\z  “  Z(x)>  c  aCj  U  az?R  u  .  (3.11) 

Since  the  more  general  situation  may  be  posed  as  a  mathematical 
problem,  we  will  say  a  few  words  about  It  at  the  end  of  this 
section.  However,  there  do  not  seem  to  be  any  Interesting  or 
Important  problems  of  flows  through  dams  that  are  excluded  by  the 
restriction  we  make  here. 

From  (3.11)  we  have 

T  =  C(x,z)  e  C«.(n)jz  =<  z*(x)>  «  u  u  a Dj  u  3nf  , 

(3.12) 

and  thus  by  using  (1.31J)  and  (3.3)  we  conclude  that  at  the  point 
(x,z*(x)),  either 

“  1  k  ’  n  »  "in'"**11  *  or  ^  *  0  .(3.13) 

The  fact  that  $  >  0  on  a for  sufficiently  smooth  3 D  (for  exam¬ 
ple,  with  bounded  curvature)  follows  by  getting  a  contradiction 
from  the  assumption  that  4>  *  0  at  any  point  of  this  set. 

Using 

r**(x) 

t(x,z)  •  4»(x ,z 1  )dr* 

*  z 

and  (3.3),  and  setting  c(x,z)  =  0  for  z  >  z  (x),  we  get 


(3.14) 


r 


AC(X,Z)  »  -  A6(x,z)  +  A6(x,Z*(x)) - ^  |jUx ,Z*(x ) ) 

n  •  k 

+  v  .(*(x,z*(x))vz*(x))  +  (3.15) 

n  •  k 

In  the  region  D  u  {(x,z)|z  >_  z*(x)}.  Here  6gD  Is  a  Dirac  measure 
on  3 D. 

On  3flf  u  (T  n  Int(3Di)) 

C  a  ?C  *  0  ,  (3.16a) 

and  just  below  this  set 

A?  =  -  — ^t|^-(x,z  (x))  ^1  ,  (3.16b) 

n  •  k 

where  we  get  equality  just  below  3flf.  Just  below  T  n  , 

(3.15)  becomes 

if 

AC  ■  1  +  V  •  (*{x,Z*(x))VZ*(x))  +  »ix».z  ,  (3.16c) 

n  .  k  30 

and  just  below  T  n  3Dj2, 

*  >  0  ,  (3.1 6d ) 

* 

AC  -  V  •  (*(x,z*(x))vz*(x))  +  ^.ZJ.X.??6  ,  (3 .1 6e) 

n  •  k  dU 

The  problem  for  c  near  as  given  by  (3.16a)  and  (3.16b), 
Is  close  to  the  steady  state  free  boundary  problem  described  In 
(1.7)  and  (1.8).  It  has  been  seen  (Ref.  3)  that  this  free  bound¬ 
ary  problem  may  be  viewed  as  the  limit  as  e  -*■  0  of  a  nonlinear 
problem  dependent  on  the  positive  parameter  e  and  defined  on  the 
whole  region  d.  Since  (3.16b)  shows  that  Ac  Is  not  necessarily 
constant  near  3n^  u  (T  n  Int(3Ci)),  the  treatment  previously 
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given  Is  modified  slightly.  Near  an^  u  (T  n  Int ) )  we  are  led 
to  consider  the  e-dependent  problem 


Ac(e)  +  m1n|— |^(x,z*(x)),-ljgc(i;(c))  -  0  ,(3.17) 


where 


9,(0  s 


C  >  e 

0  <  c  <  £ 


The  boundary  condition  Is  that  c(e)  0  as 
may  consider  the  limit  as  e  0  of  ;(e)  In 
tlon  of  the  equation 


(3.18) 

(x,z)  -*•  (x,Z(x)) .  We 
(3.17)  to  be  the  solu- 


r  ml 


\n  • 


+ 


f^(x,z  (x)),-l\ 


an 


C  >  0 

-  0  .(3.19) 

c  -  0 


By  using  monotonicity  arguments  for  elliptic  equations  we 
can  show  that.  If  c  and  c(£)  satisfy  the  same  elliptic  boundary 
value  problem  in  a  region,  with  the  same  boundary  data,  except 
for  the  replacement  of  (3.16b)  by  (3.17)  In  parts  of  the  region, 
then 


C  i  C(e)  1  C  +  e  •  (3.20) 

We  now  want  to  look  at  the  behavior  of  c(e)  where  0  ^  c(£)  £, 

and  for  x  such  that  (x,z*(x))  £  anf  u  (T  n  Int (az^ ) ) .  To  facili¬ 
tate  the  discussion,  let  us  make  the  definition,  for  a  function  u, 
of 


S£(u)  =  {(x,z)  £  £»|0  <  u  <_  £)  .  (3.21) 
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r 


It  follows  from  (3.20)  that  S£(c(e))  c  S£(c).  In  addition  we 


define 


A  =  ((x,z)  £  £>|(x,z  (x))  e  30^  u  (T  n  Int(3£>i))} 


If 


^^(x^x)) 


(3.22) 

Is  bounded  for  x  e  T  n  Int(3Z>i),  one  sees 


that  ( ac («■ )  {  bounded,  uniformly  In  e,  In  A  n  S£ (c ) ,  and  thus 
one  gets  equl continuity  of  the  derivatives  7c(e).  Hence 

vc(e)  7c  polntwlse  ,  (x,z)  c  S£(c)  n  A  .(3.23) 

Since 


sup  sup 

^  =  (-c,)  0  as  £  0  ,  (3.24) 

A  n  S£(c)  A  n  S£(c) 


(3.23)  Immediately  Implies  that 

<|>{e)  -*■  o  as  £  0  ,  (x,z)  £  S£(c(e))  n  A  .  (3.25) 

Let  the  equation  of  the  surface  c(e)  *  f  be  written  in  the 
form  z  *  Z£(x),  that  Is, 

c(£,x,Z£(x))  »  c  .  (3.26) 

Then  on  the  Intersection  of  this  surface  with  A,  It  will  follow 
from  (3.25)  and  (3.17)  that  we  have  In  the  limit  as  e  +  0 

-  ~ t  frtU.x.z  (x))  -  Ac ( e»x »z  ( x) )  - 
n  .  k  3 

-  m1n(-^^xj*-M).-l\  .(3.27) 

\n  •  k  3n  ' 

Similarly,  we  get  from  (3.25)  and  (3.17) 


^f^U.x.z)  ^  Ac (c »x ,z ) 

•  K 


(x.zWx.rrxTy 0 


(3.28) 
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Corresponding  to  (3.17)  for  c(e)  we  find  the  following 
equation  for  <K«)  near  8n^  u  (T  n  Int  (»£>i ) ) : 

rifclil  0  <_  5(e)  <  c 

,3 .29)  • 

VO  c(e)  >  e 

This  equation  Is  not  useful  to  compute  with,  except  for  certain 
special  types  of  boundary  data  and  regions  D,  such  as  those  occur¬ 
ring  in  the  simplest  of  the  problems  studied  by  Balocchl  and 
Magenes  (Ref.  2),  because  of  the  appearance,  at  points  on  dD  n  sn, 

of  the  term  f-if»(x,z  (x)),  which  Is,  after  all,  an  unknown  of  the 
9  n 

problem. 

On  the  other  hand,  difficulties  associated  with  the  fact  that 
3 D  and  an  have  a  nonempty  Intersection  should  be  no  difficulties  at 
all,  because  In  that  regard  the  free-boundary  aspect  of  the  prob¬ 
lem  disappears,  and  the  solution  of  the  problem  Is  straightforward 
when  the  region  occupied  by  the  fluid  Is  known.  Thus,  In  order  to 
see  what  freedom  we  may  have  In  using  versions  of  the  equation 
which  differ  from  (3.29)  and  may  be  more  amenable  to  numerical 
Implementation,  let  us  examine  the  case  when  3Qf  *  J). 

Suppose  we  want  to  solve 

A*  *  0  ,  X  e  D  ,  (3.30a) 

with 

i|»(x,z)  -*■  0  for  (x.z)  -►  (x,Z(x))  ,  x  c  r  ,  (3.30b) 

and  other  data  Imposed  to  make  >_  0  In  d.  This  Is  not  a  free 
boundary  problem,  so  that  (3.3f)  cannot  be  Imposed.  For  t  given 


by  (3.6), 


AC 


^pf£(x,Z(x))  ,  x  e  r  .  (3.31) 


n  •  k 


This  can  also  be  viewed  as  the  limit  as  e  -*■  0  of  the  solution  of 


ac(c)  +  — |^(x,Z(x))g^(c(e))  -  0  ,  x  £  r  , 


n  •  k  6  (3.32a) 

cU.x.z)  -►  0  for  (x.z)  -*■  (x,Z(x))  ,  x€r  , (3.32b) 

with  c(f)  >_  0  In  D.  As  before,  we  find  that  (3.25)  holds  for  all 
points  for  which  x  £  r  and  zU)  i  e.  Hence  as  £  -►  0  we  expect 


.  3vp  (e  »x  ,z  (x ) )  , 

—  AcU.x.z  (x))  ■ - 5r|^<x,z(x)) 


3n 


n  •  t  3n 


(3.33) 


In  addition.  It  follows  from  (3.25)  and  (3.32a)  that 


ac(£,x,z(x))  -  0 

n  •  k  3n 


(3.34) 


Differentiation  of  (3.32a)  with  respect  to  z  yields 


A*(e)  +  — 1 —  §*<x.Z(x)) 


0 i  c(c) < 


0  c(e)  >  e 

Consider  now  the  following  variation  on  (3.32): 


x  £  r 
(3.35) 


Ac («)  -  h(x)g  (5(e))  +  (— kp  |^<x ,Z(x) )  +  h(x))=  0  ,  x  e  r 

€  \n  •  k  3n  f  (3.36a) 

c(e»x,z)  -*■  0  for  (x.z)  s.  (x ,Z(x) )  ,  x  £  r  , (3.36b) 

with  c(€)  >_  0  In  d .  h(x)  Is  constrained  by  the  Inequalities 


0<  h(x)  <  -  — l-y|J(x,Z(x))  ,  x  £  r  .  (3.37) 

n  •  k 

The  same  argument  that  led  to  (3.25)  holds  again,  for  x  e  r  and 
c(c )  le.  Thus  we  get  as  e  -►  0 
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— ^  |^(e,x,z£(x))  -*■  A;(e,X,2£(x))  ■ 
n  •  k 

-  — ^fifr.Ztx))  .  X  C  r  ,  (3.38a) 

n  •  k  an 


- f=*(«»x,Z(x))  ♦  A5(e,x,Z(x))  = 

n  •  k  3n 


- ^  |J<x,Z(x))  -  h(x)  ,  x  e  r 

n  •  k  dn  (3.38b) 

Dlfferentatlon  of  (3.36a)  gives  us  the  following  equation  for 


♦(e): 


Ai|/(e)  -  h(x)  < 

V  0 


0  1  C(e)  < 


C(e)  >  e 


0 


x  e  r 

(3.39) 


Comparing  (3.33)  and  (3.38a),  we  observe  that  we  get  the  same 
result  In  the  limit  whether  we  use  (3.35)  or  (3.39),  as  long  as 
h(x)  satisfies  the  constraint  (3.37). 

Now  let  us  return  to  our  problem.  We  see  that  we  can  compute 
iKc),  near  30^,  u  (T  n  Int ( 3Z?i) ) ,  from  the  equation 


A<Ke) 


0  <  t(e)  <  e 

»  0 


V  0  C(e)  >  £ 


(3.40) 


For  the  more  general  situation  described  by  (3.15),  we  are  then 
led  to  consider 
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/Jllil  0  <_  c(c )  <  e 

Atp(e)  -  A«2(x,z)  -  j  ■  0  ,  (x,z)^e  • 

(  0  ZU)  >  e 

Corresponding  to  (3.41)  we  have  the  equation  for  5(c), 

At(e)  +  A0  ( X  ,Z  )  -  A6(x  ,Z(x)  )  -  g  (  c(e  )  )  +  — |r^(e  ,X  ,Z(x  )  ) 

€  n  •  k  3n 

-  V  •  (4'(«»x»Z(x))VZ(x)) - ,x,Z(x))63„  »  0  (3.42) 

n  •  k 

for  (x,z)  e  D  u  { (x ,z ) | z  >_  Z(x)} .  Here,  of  course,  c(e)  Is  given 
by 

fZ(x) 

c(c,x,z)  *  i|»(£,x,z')dz'  .  (3.43) 

Jz 

In  (3.42),  at  values  of  x  for  which  z*(x)  «  Z(x)  and  we  have 

prescribed  boundary  conditions  on  ^(x,z  (x)),  the  terms 

— kfk^x.Zlx))  -  v.(*(e,x.Z(x))vZ(x))  -  — L_^(e,x,Z(x))6 
n  •  k  30  n  •  k 

will  adjust  themselves  In  the  limit  as  e  -*-0  to  leave  the  equation, 
when  expressed  at  z  *  z£(x).  Invariant  under  a  suitable  change  In 
the  coefficient  of  g£.  However,  for  values  of  x  such  that 
Z(x)  t  z*(x),  that  is,  such  that  (x,z*(x))  e  an^,  such  an  adjust¬ 
ment  cannot  take  place,  because  In  the  limit  as  e  0  we  will 
Invariably  have 

1<(c»x,Z(x))  -*■  0  ,  — ]—+  |^(e,x,Z(x))  -*■  0  .(3.44) 

n  •  k 

Thus,  the  coefficient  of  g£  cannot  be  changed  for  such  values  of  x 
without  changing  the  problem  in  the  limit  c  -*■  0.  In  rough 
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physical  terns,  this  situation  arises  because  on  T  n  3D  we  have 
an  outflow,  under  assumption  (3.11),  and  we  can  do  whatever  we 
want  with  the  flow  as  It  leaves  the  region,  as  long  as  we  do  not 
create  disturbances  which  propagate  back  Into  the  region.  (When 
(3.11)  Is  not  satisfied  and  we  have  a  point  e  a d2  at  which 
z  (x)  ■  Z(x)  with  an  Inflow,  we  are  still  free  to  make  suitable 
adjustments  In  the  coefficient  of  g£  In  (3.42),  as  the  Inflow  Is 
effectively  prescribed  at  such  points,  through  the  boundary  condl 
tlons.  Independently  of  such  variations.)  The  limit  as  €  0  of 

t(e )  will  satisfy 

AC  +  A6(x,z)  -  A0 (x ,Z(x) )  +  — VJrj  |^|f(e,x,Z(x)) 


11m 

e+o 


(<|/(€,x,Z(x))vZ(x)) - 11m  i|»(c,x, Z(x))6 

n  •  k  €  0 


I1 

lo 


C  >  o 

-  0  ,  (3.45) 

c  »  0 


where  ^(c)  Is  given  by  (3.41)  and  appropriate  boundary  conditions 
on  3D. 

For  the  purpose  of  solving  (3.41)  and  (3.42)  numerically,  we 
Introduce  a  pseudo-time  variable  which  we  label  "t,"  but  which 
should  be  distinguished  from  the  real  time  appearing  earlier  In 
this  paper.  Thus,  we  write  In  place  of  (3.41)  and  (3.42) 

*t(e,t)  -  A*(e,t)  -  A62(x,z) 

0  <_  c(e  ,t)  <  e 
C(e,t)  >  e 


0 
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,  (x»Z,t)  £  Dx(0,«) 

(3.46) 


and 

Ct(e,t)  -  Ac(e,t)  +  A«(x,z)  -  A© (x ,Z(x) )  -  g£(cU.t) 

+  — |-^(e,t,x,z(x))  -  V  *(<|/(e  ,t,X,Z(x))?Z(x)) 
n  •  k  3n 

-  — U^(e,t,x,Z(x))6  _  ,  (x,Z,t)  6  z?x( 0,«) 

n  •  k  (3.47) 

We  are  now  ready  to  discuss  numerical  solution  of  the 
problem.  However,  before  we  proceed,  we  should  note  one  further 
complication  that  arises.  This  Is  that  In  the  problem  (3.3)  for 
<|i  In  the  region  n,  which  we  may  now  replace  by  a  problem  for  ^  1  n 
D  through  (3.4)  and  the  Interpretation  of  *!>  as  the  limit  of  <Ke) 

In  (3.41),  all  the  boundary  data  for  ^  on  3 D  are  homogeneous  with 
the  exception  of  the  data  on  30^  n  3fi,  given  by  (3.3e).  The  dif¬ 
ficulty  Is  that  we  do  not  generally  know  the  extent  of  this  set 
until  we  have  solved  the  problem.  Hence,  we  shall  assume  for  the 
immediate  future  that  30^  *  9,  and  see  how  to  solve  the  problem 
in  that  case.  Then  we  will  show  how  an  elementary  modification 
of  the  algorithm  used  for  the  solution  when  30^  =  9  will  yield 
the  solution  for  the  practically  important  case  3 +  9.  (There 
are,  however,  cases  where  we  know  a  priori  that  30^  c  3«.  in 
such  situations  we  can  make  the  appropriate  changes  In  &,  defined 
In  (3.2),  to  get,  for  ■  p  +  02,  ■  0,  x  e  30^,  Instead  of 

(3.3e).  Then  the  algorithm  which  we  give  presently  will  be  appli¬ 
cable  without  any  modification.  Such  a  case  occurs  when  30^*  BQ, 
given  by  (1 .14).) 
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As  a  preliminary  to  the  discussion  of  approximate  solutions 
of  the  dam  problem,  we  Introduce  some  operators.  In  view  of  the 
boundary  conditions  Imposed  on  $  In  (3.3)  and  our  tentative 
assumption  that  0  0,  we  start  with  the  operator  S0(t) 
defined  as  follows: 

C(t)  3  S0(t)c0 

Is  the  solution  of  the  Initial  value  problem 
Ct  a  At  .  (x.t)  e  Dx(0,«>) 

t  3  0  ,  (x.t)  e  (3 D0  u  aZ?R)x(0,») 

|^-0  ,  (x.t)  e  3Z?jX(0,«) 

?(t)  ^  Co  »  X  e  D  . 

It  follows,  then,  that  the  solution  of 

3  A*  -  A6Z  ,  (x.t)  e  Dx( 0,») 

3  0  ,  (x,t)  e  (3 Z>„  u  3^r)x(0,») 

^■*0  ,  (x.t)  e  3^jx(0,«»)  , 

uo  *0  » 

Is  given  by 

*(t)  *  S0(t)*0  -  So(t*)dt'^0z 
when  as z  Is  Independent  of  t. 

Next,  let  us  Introduce  the  operator  S(t)  defined  by 

(S(t)co )(x»z)  =  |^(X)  {(S0(t)(-50Z))(x,z' ) 

1  S0(t')df)Ae.)(x,z')jdz'  .  . 

0  /  z/  )  (3.50) 


(3.48a) 

(3.48b) 
,  (3.48c) 

(3.48d) 
(3.48e) 

,  (3.49a) 

,  (3.49b) 

(3.49c) 
(3.49d) 

(3.49e) 
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S(t)  Is  a  semi-group: 

S(tl  +  t2)  -  S(ti)S(t2)  .  (3.51) 

Remark  3.1:  In  various  problems  we  may  want  to  work  with 
one  of  the  norms 

N^u)  =  [  ju  Idxdz  (3.52a) 

>D 

or 

sup  f® 

N2(u)  s  N  1  |u_(x,z)|dz  .  (3.52b) 

x  e  R  1  Jo 

We  observe  that  the  operator  S(t)  Is  contractive  In  Nj: 

Ni(S(t)u1  -  S(t)u2)  <  Ni(Ui  -  u2)  •  (3.53) 

In  general  this  Is  not  true  in  N2.  A  case  for  which  S(t)  Is  con¬ 
tractive  with  the  norm  N2  Is  the  case  »  B0,  as  given  In  (1.14X 
Remark  3.2:  Except  for  special  cases,  such  as  that  for 

which 

(x  | ( x , z  2 )  £  Ct(p)}  3  {x  |(x ,z  i)  £  Ct(D)}  If  z2  ^  Zi  and  3Z>j  ®  B0  , 

(3.54) 

It  Is  not  generally  true  that  the  operator  S(t)  Is  monotone: 
ti(x,Z(x))  =  t2(x,Z(x))  >*0  ,  ti  >  Kz 

=4=>  (S(t)d)(x,Z(x))  -  (S(t)c2)(x,z(x))  ■  0 

S(t)ci  >_  S(t)c2 

(3.55) 

The  case  (3.54)  does  not  appear  to  be  of  practical  Interest. 
However,  we  have  a  different  type  of  monotonicity: 

Ci(x,Z(x))  ■  c2(x,Z(x))  -  0  ,  ci  >  t2  ,  cIZ  l  c2Z 

=£(S(t)c,)(x,Z(x))  -  (S(t)t2)(x,Z(x))  -  0 

s(t)c\  >  S(t)c2  .  (S(t)ci )z  <  (S(t)c2)z 

(3.56) 
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We  also  Introduce  the  operator  P(t)(t  >  0)s 


P(t);  =  c  -  t 

(3.57) 

and  the  operator  M: 

Me  =  max(e,0). 

(3.58) 

(With  this  notation,  (1.9c)  can  be  written 

vn+1  -  MP(T)S2(T)vn  . 

(3.59)) 

Remark  3.3:  It  Is  not  generally  true  that  the  operator 
MP(t)  Is  contractive  In  either  of  the  norms  N^,  1  *  1,  2,  given  In 
(3.52).  However,  If  clz  <.  0,  c2Z  <_  0,  and  ci(x,Z(x))  * 

C2 (x ,Z(x ) )  ■  0,  we  have 

N.(MP(r)Cl  -  MP(tk2)  iN.Uj  -  C2)  ,  1®1,2. 

7  7  (3.60) 

Also,  monotonicity  holds: 

Cl(x,Z(x))  *  c2  <x  ,Z(x))  -  0  ,  >  ?2  »  Clz  < 

=^(MP(t)ci)(x,Z(x))  »  (MP(t)C2)(x,Z(x))  -  0  ,  MP(t)ci 

1  MP(t)?2  ,  (MP(t)Ci)z  <  (MP(t)c2)2  . 

(3.61) 

Because  of  the  frequency  with  which  we  will  use  the 
monotonicity  relations  (3.56)  and  (3.61),  we  shall  use  a  special 
notation: 

Ci  >*  C2  and  only  If  ti(x,Z(x))  «  c2(x,Z(x))  -  0  , 

Ci  >  C2  »  and  C12  <  c22 

(3.62) 

In  words  we  will  say  that  cj  Is  more  *  than  or  equal  to  c2,  c2 
Is  less  *  than  or  equal  to  ci»  ate. 


35 


Let  us  now  turn  to  the  dam  problem,  taken  to  be  the  limit  of 
(3.46),  (3.47)  as  e  0.  We  generate  approximate  functions 
^n(x,z)  and  cn(x,z),  which  are  supposed  to  approximate  ^(x,z,nx) 
and  c(x,z,nx),  respectively: 

4*n ( x »z)  ~  *(x,z,nx)  ,  (3.63a) 

Cn(x,z)  ~  c(x,z,nx)  .  (3.63b) 

i^n  and  cn  will  be  required  to  have  the  following  properties: 

5n(x,Z(x))  -  0  ,  (3.64a) 

♦"  >  0  ,  (3.64b) 

cn  >  0  .  (3.64c) 

Our  solution  of  (3.46),  (3,47)  Is  done  by  a  "split-step" 
scheme.  Note  that  If  (3.46)  were  just 

+t  *  -  48z 

with  the  boundary  conditions  (3.3)  and  the  assumption  30  j  3  0, 
we  would  have  simply 

*((n  +  1)t)  ■  -  ^S(x)^J  ip(x  ,z‘ ,nt)dz'jjz  .  (3.65) 

In  that  case, 

C((n  +  1 )t)  *  S(x)c(nx)  (3.66) 

would  solve  (3.47)  except  for  the  term  -g£(t).  (Note  that 

c  >*  0  =>  s(t)c  >*  o  (3.67) 

as  long  as  aez  <  0,  as  required  In  (3. Zb).)  We  approximate  the 
effect  of  the  term  -g£(c).  acting  over  a  time  t,  on  c  by  the 
operator  P(x).  Then,  In  view  of  the  fact  that  we  are  to  have 
cn+1  >*  0,  we  operate  on  the  result  by  M,  since  any  violations  of 
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this  condition  can  only  be  due  to  errors  In  the  split-step  scheme. 
Thus,  our  algorithm  Is 

Cn+1  ■  MP(T)$(T)tn  -  Fn+1c0  .  (3.68a) 

where 


Fc  =  MP(t)S(t)c  •  (3.68b) 
Thus  the  procedure  for  solving  the  dam  problem  Is  just  a  varia¬ 
tion  on  algorithm  II  of  section  1  for  the  diffusion-consumption 
problem.  This  method  has  been  called  the  "truncation  method"  In 
the  literature  (Ref.  3). 

It  Is  not  hard,  given  a  particular  dam  D  and  particular 
boundary  conditions,  to  find  functions  c~  and  i \  such  that 


.  *  _ 
U  >  C 

_  +  *  + 
Ft  <  { 


(3.69a) 

(3.69b) 


(For  example,  •  0  always  works;  so  does 

.  fZ(x)  + 

5  (x,z)  ■  j  ^  (x.z')dz'  where  ♦  (x,z)  satisfies 


«  0 

3 ’ll* 
an 


A0- 


X  e  D 


X  £  3D()  U  9  Ls;-' 


x  €  a  d. 


I 


and  2(x)  may  be  taken  to  be  Z(x)  when  Z(x)  <  •;  Z  may  be  taken  as 
a  suitable  finite  bound  when  Z(x)  *  ®,  chosen  so  that  the  Integral 
for  c+  Is  finite.)  Then  the  sequence  {FnO  will  be  nondecreas¬ 


ing  *  and  the  sequence  {Fnc+}  will  be  nonincreasing  *.  Further, 

with  the  definitions 

fZ(x) 

»1n  Md.Ca)  *  »1n(-C ,,(x,z')  ,  -c  (x,z')dz* 

z  12  22  (3.70a) 
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max  *  (?ltc2) 
we  see  that 


Z(x) 


max(-c12(x,2* ) 


-4 


2Z 


(x,z'))dz'  , 

(3.70b) 


Fn  max  *  U",c+)  >*  max  *  (fY.fY)  .  (3.71) 

From  the  contractiveness  of  F  In  Nx,  as  shown  by  (3.53)  and 
(3.60),  one  may  deduce 

MfY)  <  Nx(max  *  (fY.fY))  <  Mmax  *  U'Y)) 

1  Nx(c“)  +  Nx(c+) 
n  -  (3-72) 

so  that  the  nondecreasing  *  sequence  {F  r,  }  Is  bounded  from 
above  *  and  hence  converges.  Clearly  the  nonincreasing  * 
sequence  {Fnc+}  also  converges. 

Denote  by  1  a  function  with  the  properties  "s(x,Z(x))  =  0  and 
<_  0,  and  which  Is  a  fixed  point  of  F: 

F?  «  ?  *  MP(t)S(t)c  .  (3.73) 


There  Is  at  most  one  such  c,  since  If  there  were  two,  Cx  and  c2, 
we  would  have 

Nidi  -  li)  •  Nx(MP(t)S(t)?x  -  MP(t)S(t)T2)  <  Nx(S(t)cx  -  S(t)?2). 

(3.74) 

However,  from  the  definitions  of  S(t)  and  Nx  In  (3.50)  and  (3.52a), 
respectively,  one  sees  that,  as  long  as  3 u  3£R  t  0,  there  Is  a 
constant  c(t)  >  0  such  that 

NxiS(t)ux  -  S(t)u2)  <  e"Tc(T)Nx(ux  -  u2)  .  (3.75) 

Inserting  this  Into  (3.74),  we  get  that  Nxdx  -  ?2)  ■  0,  or 
Cl  ■  5^2*  This  Is  not  an  assertion  of  uniqueness  of  the  steady 
state  solution  of  the  actual  physical  problem,  but  of  uniqueness 
of  the  approximate  steady  state  solution  generated  by  (3.68). 
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Remark  2.1  mentions  the  problem  of  uniqueness  of  the  steady  state 
solution  of  the  physical  dam  problem.  (Elementary  analysis  shows 


that,  for  some  characteristic  dimension  a  of 

c(t)  •  e'®2/T  (3.76) 

will  work.  However,  this  Is  a  very  conservative  estimate.  The 
pseudo  "time-dependent"  problem  of  (3.46)  and  (3.47),  the  solu¬ 
tion  of  which  Is  approximated  by  (3.68),  Is  different  from  the 
actual  time-dependent  problem  discussed  In  the  last  section. 
However,  as  c ( t )  In  (3.47)  begins  to  converge  to  c,  supp  c(t) 
will  vary  slowly,  and  the  operator  F  In  (3.68b)  will  look  at  late 
times  like  the  operator 

?U  2  Xsppp  c(t)S<T>“  <3-77a> 

where 

fZ(x) 

(X^Mx.z)  =  J  XE(x,z')(-u2(x,z'))dz'  (3.77b) 

N  1  ^ 

for  any  set  e  c  R  x(0,«).  Repeated  application  of  the  operator 
f  generates  a  function  c  such  that 

M-Cj)  -  A0z  ,  x  «  supp  C(t )  ,  (3.78) 

(Ref.  10),  and  thus  for  p  given  by  (3.1),  Ap  -  0  In  supp  c(t). 
Hence,  at  late  time,  the  free  boundary  does  evolve  In  a  manner 
similar  to  the  evolution  of  the  free  boundary  In  the  real  time- 
dependent  problem.  In  such  a  case,  by  remark  2.1,  we  may  expect 
an  exponential  convergence  of  cn  to  T  as  n  -*■  ”,  for  the  case  when 
all  components  of  supp  ?  are  connected  to  3 D0  u  and  we  take 
a  nondecreasing  *  sequence  to  get  to  t". ) 


i 

i 

I 


i 
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We  can  do  the  split-step  scheme  (3.68a)  generating  sn+1 
from  ?o  reverse  order  to  get 

S'n+1  =  (S(t )MP(t ))n+1£o  ,  (3.79) 

which  converges  to  s"'  as  n  It  Is  clear  that 

-  S(t)c  ,  S«MP(t)c'  ,  (3.80) 

and  thus 

?'  -  T  <  F  <  ?  .  (3.81) 


In  the  relatively  uninteresting  case  (3.54)  where  the  operator 
S(t)  Is  monotone,  we  can  show  that  the  true  solution  s  satisfies 

S  <  C  <  Z' 


and  thus  use  (3.81)  to  obtain  an  error  estimate.  This  sort  of 
estimate  was  given  already  In  (1.11)  for  algorithm  II,  and  we 
have  essentially  reproduced  Its  derivation  (Ref.  3).  Unfortu¬ 
nately,  for  the  general  case  not  covered  by  (3.54)  we  do  not  gen¬ 
erally  expect  Z'  to  be  an  upper  bound  for  S.  Indeed,  we  would 
anticipate  that  Z'  <  s  at  points  in  D  near  to  the  Interior  of 
T  n  3Dr. 

We  have  not  at  this  time  produced  an  error  estimate  for  the 
difference  between  Z  and  c,  the  solution  of  (3.45).  The  follow¬ 
ing  steps  Indicate  how  we  would  proceed  to  get  a  precise  estimate. 

Step  1.  We  observe  that  (3.73)  can  be  written,  with 
n  =  supp  c  n  D,  as 

Z  =  X.S(t)c  .  (3.82) 

n 
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Also,  we  note  that,  given  regions  and  n2  =  flj,  the  solutions 
of 

Cl  -  XniS(x)ci  ,  C2-XqS(t)c2  ,  (3.83) 

satisfy,  for  a$2  ^  0, 

* 

Cll  .  (3.84) 

Step  2.  A  check  of  (3.82)  and  (3.78)  shows  that  c”  - 
satisfies  the  correct  differential  equation  In  n.  Thus,  the  main 
task  involves  checking  to  what  extent  the  boundary  conditions  on 
e  are  satisfied  by  C.  Denoting  by  T  the  set  of  points  on  the  top 
of  an, 

_  (  _  sup  \ 

T  5  J (x ,z )  c  C*(n)|z  -  _z'(  ,  (3.85) 

(  (x  ,z 1 )  e  n  ) 

we  have  no  difficulty  In  verifying  the  condition 

?(x,z)  =•  0  ,  (x,z)  €  T  ,  (3.86) 

from  (3.73).  The  error  with  which  the  other  boundary  conditions 

hold  on  parts  of  an  near  30  may  be  ascertained  without  too  much 

difficulty.  The  only  problems  arise  on  parts  of  an  which  remain 

uniformly  removed  from  3 D  In  the  limit  t  ■+■  0.  We  label  these 

parts  of  an  as  an^ . 

Step  3.  Since  we  do  not  have  any  a  priori  estimate  of  the 
regularity  of  3n^,  we  may  consider  a  modified  algorithm  which 
always  generates  a  "smooth"  boundary.  The  deviation  of  the  solu¬ 
tion  generated  by  this  modified  algorithm  from  "c  can  be  bounded 
by  using  a  variation  on  the  monotonicity  arguments  and  operators 


Introduced  by  Brezls,  Berger,  and  Rogers  (Ref.  6)  to  study  error 
bounds  for  approximate  solutions  of  the  Stefan  problem.  We  let 

-Z(x)  sup 
'  z 


(U*u)(x,z)  2  J 


.  ,  (-u_(y))dz' 

|y  -  (x,z')|  <_ 6 


,(3.87a) 


and 


(U“u)(x,z)  =  j 


Z(x)  Inf 

(-u  (y))dz' 
!y  -  (x.z1)}  <  6  2 


(3.87b) 


For  a  set  f  c  d,  we  define 


(  ,/ 

) 

j(x,z)|^u; 

s(i2  *<*>“') ) 

(3.88) 


If  3 d  has  curvatures  bounded  In  magnitude  by  1/6,  the  set 
UjU*E  Is  the  complement  of  the  set  of  points  In  the  complement  of 
E  through  which  the  surface  of  a  closed  ball  of  radius  6  can  be 
passed  In  such  a  way  that  the  whole  ball  lies  outside  E.  The 
boundary  of  U~U*E  Is  generated  by  rolling  a  "marble"  of  radius  6 
over  the  boundary  of  E,  keeping  the  marble  outside  e. 

When  u2  0  and  u(x,Z(x))  *  0,  we  have 

MP(x)U*u  iX%W*)|u>T>“  1  "Xx<<*,0|u>,>u  S*  • 

8  -  65  ”  (3.89) 

One  may  consider.  In  place  of  (3.68),  the  modified  algorithm 

Cx"  “  Fx"?0  ’  (3.90a) 

where 

Fxc  2  XUJuJ{(x,z)1s(t)cit>S^t^  ’  (3.90b) 

n 


As  n  ■*■  *,  C  "  will  approach  a  "steady  state"  C  satisfying 

X  X 


Sc  “  XU5UjMx,z)is(T)txlT} 


r  ..\S(t)cv 


(3.91) 
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It  Is  clear  from  (3.84)  that 


_  *  _ 


V>  5  .  (3-92) 
and  moreover  ■  supp  n  D  has  the  property  that  parts  of  anx 
which  project  Into  JTx  will  be  smooth. 

To  bound  Ni(Fx  -  ?)»  we  may  formulate  a  third  algorithm, 
according  to  which  a  quantity  Ut)n  Is  constructed  as 


where 


(?*)"  a  (F?>n5o  . 


F?C  =  MP(x)(S(x)t)c 


and  the  operator  S(x)*  Is  required  to  satisfy 

S(t)X?  L*  ^S(x)c  . 


(3.93a) 


(3.93b) 


(3.93c) 


Then,  upon  setting  s0  =  0  In  (3.68),  (3.90),  and  (3.93),  we  obtain 
for  the  steady  state  c*  of  (3.93), 


-+  *  —  *  — 
^  ?x  2-  5 


(3.94) 


and  we  can  bound  Ni(?x  -  D  by  Ni(Cg  -  c"). 

We  may  find  an  appropriate  operator  S(x)*  by  replacing  A6z 
and  S0(t)  as  they  occur  In  (3.50)  and  (3.48)  by  ( A6 2 ) ^  and  So(t)*, 
respectively.  (For  example.  In  the  case  where  =  0,  we  may  let 
D  ■*  U*d  In  the  definition  of  So(t)*  and  we  may  let  (Ae)£  ■  U^Ae.) 

+  i  .+  * 

Given  S0(x)fi  and  (Aez)fi,  one  can  use,  for  u2  >.  0, 

MP(x)(ui  +  u2)  <*  MP(x)ui  +  u2  ,  (3.95) 

to  get 
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(Fj)n+1(0)  -  Fn+1(0)  <*  j^X>  j^+l)T  [SO(t)J(-A0z)J 

-  S0(t)(-A6  )]dtdz' 

*  f Z(x )  r“  .  . 

<  ]z  Jo  cs0(t);(-.ez); 

-  S0 (t ) (-A0  } Idtdz  *  . 

Z_+  _  (3-96) 

This  last  bound  Is  seen  to  be  an  upper  *  bound  on  ;  •  ?  >*  0. 

6  — 

It  may  be  evaluated  directly  by  observing  from  (3.48)  that 

T  »  f  S0 { t )d tip 
■*o 

satisfies  (if  3 d0  u  aZ?R  /  0) 

M  ,  X  e  D  ,  (3.97a) 

V  =  0  ,  x  £  3D0  u  3Z)r  ,  (3.97b) 

f^=0  ,  X  e  30j  .  (3.97c) 

Finally,  one  may  obtain 

Mi*  -  ? )  <.  C-  (lg)z  +  ?z]dxdz  ,  (3.98) 

where  D  Is  a  suitable  region  for  which  the  Integral  can  be 

bounded  and  which  we  know  a  priori  to  contain  supp  i*. 

Step  4.  It  Is  obvious  from  (3.91)  that  "ix  =  0  on  the  set 

i  sup  \ 

T  *  (x,r)  e  C*(n  )Jz  •  z*  .  (3.99) 

x  *  x  (x,z')  e  nx  ( 

It  can  also  be  shown  (Ref.  10)  that  for  given  6,  and  x  -*•  0, 

(cx)z  ■  o((x  in  if)  on  {(x,z)  «  Tx  |d1st((x,z) ,30R  u  a >  6> 

(3.100) 

Regarding  the  free  boundary  condition,  that  4^  -*■  1  just  below 
the  set  ((x,z)  e  Tx ]d1st((x,z),3Z?)  _>  6),  that  will  follow,  In  the 
limit  as  x  0,  directly  from  (3.91). 
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Our  guess  ts  that  careful  handling  of  all  the  estimates  will 
finally  yield,  for  the  Ni  norm  of  the  difference  between  ?  and  c, 
a  result  which  Is  o|t  in 

To  this  point  we  have  assumed  that  30^  *  0.  Let  us  now  see 
how  we  may  modify  the  algorithm  (3.68)  when  30^  +  0.  The  follow 
Ing  scheme  may  be  used. 

$o(t)  Is  still  given  by  (3.48).  Let  Go(t,x,z,x,z)  be  the 
derivative  of  ($o(t)Co)(x,z)  with  respect  to  Eo(x,z): 

(S0(t)C0)(x,z)  »  I  Go(t,x,z,x,z)£o(x,z)dxdz  •  (3.101) 

We  next  Introduce,  in  place  of  S(t)  In  (3.50),  the  operator 

(S*(tko)(x,z)  =  £(X)  |(S0(t)(-C02))(x.z*) 


-((j  So(t')df)aez)(x,z') 

-  |  (|  G0(t'  ,x,z'  ,x,z)dt'jk  •  n(x,z)ds|dz' 

1  1  li  inn 


(3.102) 

S  (t)  Is  the  operator  that  would  be  appropriate  In  the  computa¬ 
tional  algorithm  If  30  =  .  In  general  this  Is  not  the  case, 

ft  ft  ft 

and  It  does  not  generally  follow  from  D  that  S  (t)co  >_  0. 
Since  we  want  the  Iterative  solutions  we  generate  to  satisfy 
tn  >  0  (nonnegative  pressure),  we  shall  replace  (3.68)  In  the 


case  3  D..  f  0  by 


where 


-*nc 
F  to 


(3.103a) 


(3,103b) 


F*C  s  MP(t)mV(t)c 
and 

M%  =  max  *  (c.0)  .  (3.103c) 

Remark  3.4:  The  various  algorithms  (3.68),  (3.90),  (3.93), 
and  (3.103)  may  be  replaced  by  algorithms  where  In  place  of 
S0 ( t )  given  by  (3.48),  we  work  solely  with  semi-groups  for  the 

Id  % 

diffusion  equation  In  R  x(0,«).  Operations  with  such  semi¬ 
groups  are  Interspersed  with  multiplication  by  operators  X„  for 
given  sets  E,  defined  In  (3.77b),  to  approximately  satisfy  the 
Dlrlchlet  conditions  on  3 D,  and  also  with  multiplication  by  suit¬ 
able  reflection  operators,  to  approximately  satisfy  the  Neumann 
conditions  on  the  remainder  of  3 D. 

For  the  sake  of  mathematical  completeness,  we  will  discuss 
the  case  where  the  set  {z^(x)l  defined  in  (3.8)  may  contain  more 
than  one  element  for  a  given  x,  and  where  the  condition  (3.11) 
need  not  hold.  The  essential  difference  In  our  approach  will 
Involve  replacing  the  operator  P(t)  In  (3.103)  by  something  more 
complicated. 

Given  quantities  Zi(*).  z2(»).  and  a  function  u(»,z),  we 
define,  for  z^*)  >  z2(*)  ^  0, 

10 

z  »  Zj(*) 

m1n(u(*,z)  -  u(*,Zx(*)),t)  .  (3.104) 

*i(‘)  tz  *  z2<‘) 

m1n(u(*,z2(*))  -  u(*,zi(*)).t) 
z2(  • )  1  z  >  0 
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We  recall  that  we  defined  the  sets  of  Integers  S^(x)  In  (3.10a) 
and  (3.10b),  according  to  whether  n(x)  was  odd  or  even. 

Let  the  operator  Q(t,x)  be  defined  by 
(Q(t,x)u)(x,z)  =  u(x,z)  -  v(-,z1(x),z;u(x,z)) 


u(z2l(x),z2l+1(x),z;u(x,z)) 


1«SD(x) 


(3.105a) 


for  Z(x)  *  •,  and  by 


(Q(t,x)u)(x,z)  =  -  ^z^.^xJ.z^txJ.zjutx.z))  +  u(x,z) 

n?0(«)  (3-105b) 

when  Z(x)  <  «.  Finally,  In  the  general  case  we  replace  (3.103)  by 

CQn  -  FQnCo  ,  (3.106a) 

where 

Fqc  h  Q(t,x)hV(t)c  •  (3.106b) 


4.  A  Different  Problem 


We  shall  consider  a  variation  on  the  problems  discussed 

heretofore,  but  only  for  the  case  3D2  '  1,  where  3D2  Is  defined  by 

(1.20)  and  (1.21).  Suppose  the  boundary  values  of  p,  prescribed 

as  pR  on  3D  In  accordance  with  (1.13d),  are  Increased  to 

P  ■  Pr*(3Dr)  +  P'  »  (x»z)  c  3°o  u  3°R  t  (4.1a) 

p'  >  0  ,  { (x ,z )  e  9D|p«  >  0)  5  r  c  3P0  u  3D  . 

R  (4.1b) 

If  we  restrict  our  attention  to  components  of  fl  such  that 
an  n  (3D0  u  3dr)  f  0,  we  obtain 

n  =>  n  ,  p'  s  p  -  p  >  0  on  3tif  .  (4.2) 
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On  the  other  hand, 

p'  ■  0  ,  (x,z)  £  (3Z)o  u  30R)  n  38  -  T  f  (4.3a) 

|{p-  0  ,  (X.Z)  e  3 Dx  n  3n  .  (4.3b) 

We  get  p'  >  0  on  3Z?j  n  an,  since  otherwise  there  would  be  a  viola¬ 
tion  of  (4.3b),  and  thus  from  (4.1)-(4.3), 

0  .  (4.4) 


Df^lS  < 

(aD«u3DB)n3n-r  3 


(3Z)0u  3Dp)n3fl- 

Accordlng  to  (1.31a),  $  ■  -  p  -  z,  ?  *  -  p  -  z,  and 

4>*  «  ?  -  =  -  p'  ,  x  *  8  .  (4.5) 


Thus 


|^-dS  >  0  . 


(4.6) 


(4.7) 


■f(3D0u3OR)n3n-r  3n 

Since  -  0  for  (x  ,z)  e  3Dj  n  3ft,  (4.6)  Implies 

f  1  0  • 

J3Dn38-T  3 

Now  we  use  the  facts  that  4$'  ■  0  in  n  and  ■  0  on  38^.: 

f  !$-«  *  I  i  °  •  (4-8> 

Jrnan  3  J3nf  3 

*  »v 

In  ij  -  n  we  have  4$  =  0  and  the  boundary  condition  |£  *  0  on 
38f  leads  to 


(4.9) 


I  |i<ls  “  f  ~  ||dS  . 

•*3nf  3n  Ja(n-n)na d  an 

Because  of  the  condition  (1.31J),  the  assumption  3 D2  •  9,  and  the 

fact  that  ■  0  on  a(S  -  n)  n  3Dj,  we  get 

[  ,  tins  ».  0  .  (4.10) 

Mfl-n)nBD-r  9n 
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Thus  one  derives  from  (4.8)-(4.10)  that 

f  !*AiS  +  f  f±dS<0  .  (4.11) 

Jrnan  3n  Ja(n-n)nr  30 

In  physical  language,  (4.11)  establishes  the  monotone 
increasing  dependence  of  the  flow  Into  the  dam  across  r  on  the 
pressures  prescribed  on  r.  Now  we  may  Imagine  a  problem  In  which 

r  1$  the  boundary  of  D  with  a  reservoir,  whose  height  is  not 

known  but  across  which  a  total  flux  Is  prescribed.  Part  of  the 
problem  Is  to  determine  the  height  of  the  reservoir.  If  one 
establishes  upper  and  lower  limits  for  the  reservoir  height,  such 
that  at  the  upper  limit  the  flow  out  of  the  reservoir  will  be  too 
large  and  at  the  lower  limit  the  flow  out  will  be  too  small,  then 
the  monotonicity  just  deduced  establishes  the  existence  of  a 
unique  reservoir  height  for  which  the  flux  condition  Is  satisfied. 
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